# executed by mmendozaavina@fas.harvard.edu on 2025-12-16
# see environment information pasted at the bottom


library(tidyverse)
library(estimatr)
library(broom.helpers)
library(modelsummary)
library(marginaleffects)
library(kableExtra)
library(ggthemes)
library(patchwork)
library(ggh4x)

# import datasets ----------------------------------------------------------------------------------------
d <- read_rds("appb.rds")
d1 <- read_rds("study1.rds")
d2 <- read_rds("study2.rds") 
d3 <- read_rds("study3.rds")


# Appendix B ---------------------------------------------------------------------------------------------

# Figure B.1
p1 <- d %>%
  pivot_longer(starts_with("sentiment_"), names_to = "dv", values_to = "y") %>%
  mutate(dv = sub("^sentiment_", "", dv)) %>%
  group_by(pid, dv) %>%
  group_modify(~tidy(lm_robust(y ~ 1, weights = weight, data = .))) %>%
  ggplot(aes(dv, estimate,
             color = pid, shape = pid)) +
  geom_point(size = 2) +
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error,
                    ymax = estimate + 1.96 * std.error),
                width = 0) +
  geom_errorbar(aes(ymin = estimate - 1.645 * std.error,
                    ymax = estimate + 1.645 * std.error),
                width = 0, linewidth = 1) +
  scale_y_continuous(limits = 0:1, breaks = seq(0, 1, .2)) +
  labs(y = "Proportion of responses",
       x = "Sentiment",
       color = NULL, shape = NULL) +
  scale_color_grey() +
  scale_shape_manual(values = c("triangle", "square")) +
  ggtitle("A. Sentiment") +
  theme_tufte(base_size = 12) +
  theme(legend.position = "bottom",
        legend.direction = "horizontal",
        legend.box.margin = margin(-10, -10, -10, -10),
        legend.spacing.y = unit(0, "line"))

p2 <- d %>%
  pivot_longer(starts_with("misinfo_"), names_to = "dv", values_to = "y") %>%
  mutate(dv = sub("^misinfo_", "", dv)) %>%
  group_by(pid, dv) %>%
  group_modify(~tidy(lm_robust(y ~ 1, weights = weight, data = .))) %>%
  ggplot(aes(dv, estimate,
             color = pid, shape = pid)) +
  geom_point(size = 2) +
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error,
                    ymax = estimate + 1.96 * std.error),
                width = 0) +
  geom_errorbar(aes(ymin = estimate - 1.645 * std.error,
                    ymax = estimate + 1.645 * std.error),
                width = 0, linewidth = 1) +
  scale_y_continuous(limits = 0:1, breaks = seq(0, 1, .2)) +
  labs(y = "Proportion of responses", 
       x = "Misinformation",
       color = NULL, shape = NULL) +
  scale_color_grey() +
  scale_shape_manual(values = c("triangle", "square")) +
  ggtitle("B. Misinformation") +
  theme_tufte(base_size = 12) +
  theme(legend.position = "bottom",
        legend.direction = "horizontal",
        legend.box.margin = margin(-10, -10, -10, -10),
        legend.spacing.y = unit(0, "line"))


p1 + p2 + plot_layout(widths = c(3, 2)) + plot_layout(guides = "collect") & theme(legend.position = 'bottom')



# study 1 Descriptives -----------------------------------------------------------------------------------

# Table 1 
facts <- d1 %>%
  filter(immigration_t == 1) %>%
  mutate(t1_q1 = ifelse(t1_q1 == FALSE, 1, 0),
         t1_q2 = ifelse(t1_q2 == TRUE, 1, 0),
         t1_q3 = ifelse(t1_q3 == FALSE, 1, 0),
         t1_q4 = ifelse(t1_q4 == FALSE, 1, 0),
         t1_q5 = ifelse(t1_q5 == TRUE, 1, 0))

rbind(facts %>%
        drop_na(bi_pid) %>%
        group_by(bi_pid) %>%
        summarize(across(starts_with("t1"), 
                         ~weighted.mean(.x, w = weight) * 100)),
      facts %>%
        summarize(across(starts_with("t1"),
                         ~weighted.mean(.x, w = weight) * 100)) %>%
        mutate(bi_pid = "Overall")) %>%
  pivot_longer(cols = starts_with("t1_"),
               names_to = "question") %>%
  mutate(value = round(value),
         claim = case_when(
           question == "t1_q1" ~ "Throughout the pandemic, most illegal immigrants received COVID-19 stimulus checks from the
federal government. [FALSE]",
question == "t1_q2" ~ "Illegal immigrants are not eligible for federal assistance programs such as Social Security, Medicare,
and Unemployment Insurance even though they financially contribute to them. [TRUE]",
question == "t1_q3" ~ "Under the Affordable Care Act, illegal immigrants can enroll in health insurance programs. [FALSE]",
question == "t1_q4" ~ "Illegal immigrants do not pay taxes because they do not have Social Security Numbers. [FALSE]",
question == "t1_q5" ~ "In the US, there is no law barring illegal immigrants from starting their own businesses and hiring
employees. [TRUE]")
  ) %>%
  pivot_wider(names_from = "bi_pid",
              values_from = "value") %>%
  arrange(Overall) %>%
  select(-question,
         Statement = claim) 

# Table 2
narratives <- d1 %>%
  filter(immigration_t == 2) %>%
  select(starts_with("t2_"),
         bi_pid, weight)

rbind(narratives %>%
        drop_na(bi_pid) %>%
        group_by(bi_pid) %>%
        summarize(across(starts_with("t2"), 
                         ~weighted.mean(.x == "WIN", w = weight) * 100)),
      narratives %>%
        select(-bi_pid) %>%
        summarize(across(starts_with("t2"), 
                         ~weighted.mean(.x == "WIN", w = weight) * 100)) %>%
        mutate(bi_pid = "Overall")) %>%
  rename("The day after Mia arrived to the US as a young adult, she began a job as a factory line worker. She regularly works over time, often on night shifts. ``I work like a horse,'' she says, ``but I am so grateful. I want the work.''" = t2_q1,
         "Nadia is a trained surgeon who left her native country with her 3 children. Over the last 10 years, she’s been an in-home nurse for a caregiving agency. ``I gave up my dreams so my children could grow up in America. And I would do it again in a heartbeat.''" = t2_q2,
         "When Victor arrived to the US, he began working as a street food vendor. Today, 25 years later, he owns his own business—a bakery—which employs 12 other people." = t2_q3,
         "In his native country, Julian was a software engineer. Now, in the US, he works full-time in a restaurant kitchen. Every day after work, he tries to spend some time on a programming sideproject. He dreams of launching a tech startup." = t2_q4,
         "Ana was a schoolteacher in her home country. Today, she works as a maid for a cleaning agency. ``I work so much that my body always hurts. And proudly so. I provide for myself and my family.'' If everything goes well, Ana plans to start her own cleaning agency in the coming months." = t2_q5) %>%
  pivot_longer(cols = !bi_pid) %>%
  mutate(value = round(value)) %>%
  pivot_wider(names_from = "bi_pid",
              values_from = "value") %>%
  arrange(Overall) %>%
  rename(Story = name)


# fact-check scores
d1 %>% summarize(mean(score, na.rm = T))

d1 %>% 
  group_by(bi_pid) %>% 
  summarize(mean(score, na.rm = T))

# outcome distributions (Table A.1)
wmean_dv <- function(x){
  x %>%
    summarize(across(starts_with("DV"), 
                     ~round(weighted.mean(.x, w = weight), 2)))
}

rbind(
  d1 %>%
    drop_na(bi_pid) %>%
    group_by(bi_pid) %>%
    wmean_dv(),
  d1 %>%
    wmean_dv() %>%
    mutate(bi_pid = "Overall")
) %>%
  pivot_longer(cols = starts_with("dv"), 
               names_to = "DV") %>%
  mutate(DV = case_when(DV == "dv_contribute" ~ "Illegal immigrants contribute",
                        DV == "dv_legalize" ~ "Path to citizenship",
                        DV == "dv_levels" ~ "Increase legal immigration")) %>%
  pivot_wider(names_from = bi_pid,
              values_from = value)

# correlations between outcomes (Table A.2)
d1 %>% 
  select(starts_with("dv_")) %>%
  select(`Illegal immigrants contribute` = dv_contribute,
         `Path to citizenship` = dv_legalize,
         `Increase legal immigration` = dv_levels) %>%
  datasummary_correlation(method = "pearspear",
                          output = "data.frame")


# Distribution of fact-checking scores (Figure A.1)
scores <- d1 %>%
  filter(immigration_t == 1) %>%
  select(score = score, 
         bi_pid,
         weight) %>%
  mutate(`0` = ifelse(score == 0, 1, 0),
         `1` = ifelse(score == 1, 1, 0),
         `2`= ifelse(score == 2, 1, 0),
         `3` = ifelse(score == 3, 1, 0),
         `4` = ifelse(score == 4, 1, 0),
         `5` = ifelse(score == 5, 1, 0)) 

rbind(scores %>%
        drop_na(bi_pid) %>%
        group_by(bi_pid) %>%
        summarize(across(!weight,
                         ~weighted.mean(.x, w = weight) * 100)),
      scores %>%
        select(-bi_pid) %>%
        summarize(across(!weight,
                         ~weighted.mean(.x, w = weight) * 100)) %>%
        mutate(bi_pid = "Overall")) %>%
  select(-score) %>%
  pivot_longer(cols = !bi_pid) %>%
  mutate(bi_pid = factor(bi_pid,
                         levels = c("Democrats",
                                    "Republicans",
                                    "Overall"))) %>%
  ggplot(aes(name, value, 
             fill = bi_pid)) +
  geom_bar(stat = "identity") +
  facet_wrap(bi_pid ~ .) +
  geom_text(aes(label = gsub(" ", "", paste(round(value),"%"))),
            vjust = -.5,
            size = 2.5) +
  scale_y_continuous(name = "",
                     limits = c(0, 40)) +
  scale_x_discrete(name = "5-point fact-checking score\n(# correct guesses)") +
  scale_fill_grey() +
  theme_minimal() +
  theme(text = element_text(family = "serif"),
        panel.grid = element_blank(),
        strip.text = element_text(size = 11),
        axis.text = element_text(size = 11),
        legend.position = "none") 



# Study 1 results ----------------------------------------------------------------------------------------

d1_long <- pivot_longer(d1, cols = c(dv_contribute, dv_legalize, dv_levels), names_to = "dv")

results_study1 <- function(subset) {
  group1 <- levels(factor(d1 %>% select({{subset}}) %>% pull()))[1]
  group2 <- levels(factor(d1 %>% select({{subset}}) %>% pull()))[2]
  subset1 <- d1 %>% filter({{subset}} == group1)
  subset2 <- d1 %>% filter({{subset}} == group2)
  
  if (!is.na(group2)) {
    mod_contribute <- list(
      lm_robust(dv_contribute ~ immigration_t, weights = weight, data = d1),
      lm_robust(dv_contribute ~ immigration_t, weights = weight, data = subset1),
      lm_robust(dv_contribute ~ immigration_t, weights = weight, data = subset2))
    mod_legalize <- list(
      lm_robust(dv_legalize ~ immigration_t, weights = weight, data = d1),
      lm_robust(dv_legalize ~ immigration_t, weights = weight, data = subset1),
      lm_robust(dv_legalize ~ immigration_t, weights = weight, data = subset2))
    mod_levels <- list(
      lm_robust(dv_levels ~ immigration_t, weights = weight, data = d1),
      lm_robust(dv_levels ~ immigration_t, weights = weight, data = subset1),
      lm_robust(dv_levels ~ immigration_t, weights = weight, data = subset2))
    names(mod_contribute) <- c("Overall", group1, group2)
    names(mod_legalize) <- c("Overall", group1, group2)
    names(mod_levels) <- c("Overall", group1, group2)
  } else {
    mod_contribute <- list(
      lm_robust(dv_contribute ~ immigration_t, weights = weight, data = d1),
      lm_robust(dv_contribute ~ immigration_t, weights = weight, data = subset1))
    mod_legalize <- list(
      lm_robust(dv_legalize ~ immigration_t, weights = weight, data = d1),
      lm_robust(dv_legalize ~ immigration_t, weights = weight, data = subset1))
    mod_levels <- list(
      lm_robust(dv_levels ~ immigration_t, weights = weight, data = d1),
      lm_robust(dv_levels ~ immigration_t, weights = weight, data = subset1))
    names(mod_contribute) <- c("Overall", group1)
    names(mod_legalize) <- c("Overall", group1)   
    names(mod_levels) <- c("Overall", group1)
  }
  
  bi <- d1_long %>%
    drop_na({{subset}}) %>%
    group_by({{subset}}, dv) %>%
    group_modify(~ lm_robust(value ~ immigration_t,
                             weights = weight,
                             data = .) %>%
                   tidy_and_attach() %>%
                   tidy_remove_intercept()) %>%
    rename(bi = {{subset}})
  
  overall <- d1_long %>%
    group_by(dv) %>%
    group_modify(~ lm_robust(value ~ immigration_t,
                             weights = weight,
                             data = .) %>%
                   tidy_and_attach() %>%
                   tidy_remove_intercept()) %>%
    mutate(bi = "Overall")
  
  fig <- rbind(bi, overall) %>%
    mutate(dv = case_when(dv == "dv_contribute" ~ "Illegal immigrants contribute \n(control mean: 0.57)",
                          dv == "dv_legalize" ~ "Path to citizenship \n(control mean: 0.67)",
                          dv == "dv_levels" ~ "Increase legal immigration \n(control mean: 0.49)"),
           dv = factor(dv, levels = c("Illegal immigrants contribute \n(control mean: 0.57)",
                                      "Path to citizenship \n(control mean: 0.67)", 
                                      "Increase legal immigration \n(control mean: 0.49)")),
           term = ifelse(term == "immigration_t2", "Narratives", "Facts"),
           term = fct_relevel(term, "Narratives"),
           bi = factor(bi),
           bi = fct_relevel(bi, "Overall")) %>%
    ggplot(aes(estimate, term,
               color = bi,
               shape = bi,
               xmin = conf.low,
               xmax = conf.high)) +
    geom_vline(xintercept = 0,
               linetype = "dotted") +
    geom_point(size = 2, 
               position = position_dodge(-.5)) +
    geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                      xmax = estimate + 1.96 * std.error),
                  position = position_dodge(-.5),
                  width = 0) +
    geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                      xmax = estimate + 1.645 * std.error),
                  position = position_dodge(-.5),
                  width = 0, linewidth = 1) +
    facet_grid(. ~ dv) +
    scale_x_continuous(breaks = seq(0, .2, .1)) +
    labs(x = "Treatment effect estimate with 90% and 95% confidence intervals",
         y = NULL, color = NULL, shape = NULL) +
    scale_color_grey() +
    theme_tufte(base_size = 12) +
    theme(legend.position = "bottom",
          legend.direction = "horizontal",
          legend.box.margin = margin(-10, -10, -10, -10),
          legend.spacing.y = unit(0, "line"))
  
  tab <- modelsummary(
    list("DV: Illegal immigrants contribute" = mod_contribute,
         "DV: Path to citizenship" = mod_legalize,
         "DV: Increase legal immigration" = mod_levels),
    coef_map = c("(Intercept)" = "(Intercept)",
                 "immigration_t1" = "Facts",
                 "immigration_t2" = "Narratives"),
    estimate = "{estimate} ({std.error})",
    statistic = NULL,
    shape = "rcollapse",
    gof_map = c("nobs"),
    output = "markdown")
  
  return(list(fig, tab))
}

# Main results (Figure 1 and Table A.3)
results_study1(bi_pid)

# Results for Independents (Figure A.2)
results_study1(independents)

# Results by education level (Figure A.3)
results_study1(bi_educ)

# Results by household income (Figure A.4)
results_study1(bi_income)

# Results by race (Figure A.5)
results_study1(bi_race)

# Results by nationality (Figure A.6)
results_study1(bi_immstat)

# Results by ideology (Figure A.7)
results_study1(bi_ideo)


# Coefficient equality tests (Table A.3)

mod_contribute <- list(
  "Overall" = lm_robust(dv_contribute ~ immigration_t, weights = weight, data = d1),
  "Democrats" = lm_robust(dv_contribute ~ immigration_t, weights = weight, data = subset(d1, bi_pid == "Democrats")),
  "Republicans" = lm_robust(dv_contribute ~ immigration_t, weights = weight, data = subset(d1, bi_pid == "Republicans"))
)

mod_legalize <- list(
  "Overall" = lm_robust(dv_legalize ~ immigration_t, weights = weight, data = d1),
  "Democrats" = lm_robust(dv_legalize ~ immigration_t, weights = weight, data = subset(d1, bi_pid == "Democrats")),
  "Republicans" = lm_robust(dv_legalize ~ immigration_t, weights = weight, data = subset(d1, bi_pid == "Republicans"))
)

mod_levels <- list(
  "Overall" = lm_robust(dv_levels ~ immigration_t, weights = weight, data = d1),
  "Democrats" = lm_robust(dv_levels ~ immigration_t, weights = weight, data = subset(d1, bi_pid == "Democrats")),
  "Republicans" = lm_robust(dv_levels ~ immigration_t, weights = weight, data = subset(d1, bi_pid == "Republicans"))
)


hyp_test <- function(mod, dv){
  map2(mod, names(mod), function(x, y){
    coef <- tibble(x$coefficients)
    t1 <- pull(coef[2,])
    t2 <- pull(coef[3,])
    
    hypotheses(x, hypothesis = "immigration_t2 = immigration_t1") %>%
      tidy_broom() %>%
      mutate(Model = y, DV = dv,
             Facts = t1,
             Narratives = t2)
  }) %>% bind_rows() %>%
    select(DV, Model, Facts, Narratives, Difference = estimate, t = statistic, p = p.value)
}

map2(list(mod_contribute, mod_legalize, mod_levels), 
     c("Illegal immigrants contribute", "Path to citizenship", "Increase legal immigration"),
     hyp_test) %>%
  bind_rows() %>%
  kbl(format = "markdown", digits = 3)


# Composite index as the outcome 

mod_index <- list(
  "Overall" = lm_robust(dv_index ~ immigration_t, weights = weight, data = d1),
  "Democrats" = lm_robust(dv_index ~ immigration_t, weights = weight, data = subset(d1, bi_pid == "Democrats")),
  "Republicans" = lm_robust(dv_index ~ immigration_t, weights = weight, data = subset(d1, bi_pid == "Republicans"))
)

## Table A.5
modelsummary(
  mod_index,
  coef_map = c("(Intercept)" = "(Intercept)",
               "immigration_t1" = "Facts",
               "immigration_t2" = "Narratives"),
  stars = TRUE,
  gof_map = c("nobs"),
  escape = FALSE,
  output = "markdown")

## Table A.6
map2(list(mod_index), 
     c("Composite index (0-1)"),
     hyp_test) %>%
  bind_rows() %>%
  kbl(format = "markdown", digits = 3)




# Study 2 results ----------------------------------------------------------------------------------------

clean_conjoint <- function(x) {
  x %>% 
    mutate(term = str_remove(term, variable),
           variable = fct_relevel(variable, "status", "country", "skill", "job", "taxes"),
           term = coalesce(term, variable),
           term = fct_relevel(term,
                              "status", "Naturalized citizen", "Green card holder", "Illegal alien",
                              "country", "Ireland", "Poland", "South Korea", "India", "Mexico", "Iraq",
                              "skill", "No high school", "High school graduate", "College degree", "Graduate degree",
                              "job", "Currently unemployed", "Works 20 hours per week", "Works 40 hours per week", "Works 60 hours per week", 
                              "taxes", "$500", "$7,500", "$25,000"))
}

overall_tidy <- lm_robust(selected ~ status + country + skill + job + taxes,
                          weights = weight, clusters = id, data = d2) %>%
  tidy_and_attach() %>%
  tidy_add_reference_rows() %>%
  tidy_add_estimate_to_reference_rows() %>%
  tidy_remove_intercept() %>% 
  tidy_add_header_rows() %>%
  clean_conjoint() %>%
  select(term, variable, estimate, std.error, conf.low, conf.high) %>%
  mutate(bi = "Overall")

results_study2 <- function(subset, status_subset = FALSE){
  
  if (status_subset) {
    eq <- "selected ~ country + skill + job + taxes"
  } else {
    eq <- "selected ~ status + country + skill + job + taxes"
  }
  
  subsets <- d2 %>%
    drop_na({{subset}}) %>%
    group_by({{subset}}) %>%
    group_modify(~ lm_robust(as.formula(paste(eq)),
                             clusters = id, weights = weight, data = .) %>%
                   tidy_and_attach() %>%
                   tidy_add_reference_rows() %>%
                   tidy_add_estimate_to_reference_rows() %>%
                   tidy_remove_intercept() %>%
                   tidy_add_header_rows())  %>% 
    clean_conjoint() %>%
    select(term, variable, estimate, std.error, conf.low, conf.high, bi = {{subset}}) 
  
  fig <- rbind(subsets, overall_tidy) %>%
    mutate(bi = factor(bi),
           bi = fct_relevel(bi, "Overall"),
           variable = fct_relevel(variable, "status")) %>%
    ggplot(aes(estimate, reorder(term, desc(term)),
               xmin = conf.low, xmax = conf.high,
               color = bi, shape = bi)) +
    geom_vline(xintercept = 0,
               linetype = "dotted") +
    geom_point(size = 2, 
               position = position_dodge(-.75)) +
    geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                      xmax = estimate + 1.96 * std.error),
                  position = position_dodge(-.75),
                  width = 0) +
    geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                      xmax = estimate + 1.645 * std.error),
                  position = position_dodge(-.75),
                  width = 0, linewidth = 1) +
    facet_grid(variable ~ ., 
               scales = "free_y", space = "free") +
    scale_y_discrete(labels = c("status" = expression(bold(`Immigration status`)),
                                "country" = expression(bold(`Country of origin`)),
                                "skill" = expression(bold(`Education level`)),
                                "job" = expression(bold(`Current employment`)),
                                "taxes" = expression(bold(`Total taxes paid in 2021`)))) +
    scale_x_continuous(name = "Effect on Pr(Choose Immigrant) with 90% and 95% confidence intervals",
                       breaks = seq(-.5, .5, .1)) +
    labs(y = NULL, color = NULL, shape = NULL) +
    scale_color_grey() +
    theme_tufte(base_size = 12) +
    theme(legend.position = "bottom",
          legend.direction = "horizontal",
          legend.box.margin = margin(-10, -10, -10, -10),
          legend.spacing.y = unit(0, "line")) +
    theme(strip.background = element_blank(),
          strip.text = element_blank())
  
  group1 <- levels(factor(d2 %>% select({{subset}}) %>% pull()))[1]
  group2 <- levels(factor(d2 %>% select({{subset}}) %>% pull()))[2]
  subset1 <- d2 %>% filter({{subset}} == group1)
  subset2 <- d2 %>% filter({{subset}} == group2)
  
  if (is.na(group2)) {
    mod <- list(
      m1 = lm_robust(selected ~ factor(status) + factor(country) + factor(skill) + factor(job) + factor(taxes),
                     clusters = id, weights = weight, data = d2),
      m2 = lm_robust(selected ~ factor(status) + factor(country) + factor(skill) + factor(job) + factor(taxes),
                     clusters = id, weights = weight, data = subset1)
    )
  } else if (status_subset) {
    mod <- list(
      m1 = lm_robust(selected ~ factor(status) + factor(country) + factor(skill) + factor(job) + factor(taxes),
                     clusters = id, weights = weight, data = d2), 
      m2 = lm_robust(selected ~ factor(country) + factor(skill) + factor(job) + factor(taxes),
                     clusters = id, weights = weight, data = subset1),
      m3 = lm_robust(selected ~ factor(country) + factor(skill) + factor(job) + factor(taxes),
                     clusters = id, weights = weight, data = subset2)
    )
    
  } else {
    mod <- list(
      m1 = lm_robust(selected ~ factor(status) + factor(country) + factor(skill) + factor(job) + factor(taxes),
                     clusters = id, weights = weight, data = d2),
      m2 = lm_robust(selected ~ factor(status) + factor(country) + factor(skill) + factor(job) + factor(taxes),
                     clusters = id, weights = weight, data = subset1),
      m3 = lm_robust(selected ~ factor(status) + factor(country) + factor(skill) + factor(job) + factor(taxes),
                     clusters = id, weights = weight, data = subset2)
    )
  }
  
  if (is.na(group2)) {
    names(mod) <- c("Overall", group1)
    
    rows <- tribble(~term, ~m1, ~m2,
                    "Immigration Status", "", "",
                    "Naturalized citizen", "-", "-",
                    "Country of Origin", "", "",
                    "Ireland", "-", "-",
                    "Education Level", "", "",
                    "No high school", "-", "-",
                    "Current Employment", "", "",
                    "Currently unemployed", "-", "-",
                    "Total Taxes Paid in 2021", "", "",
                    "$500", "-", "-")
  } else {
    names(mod) <- c("Overall", group1, group2)
    
    rows <- tribble(~term, ~m1, ~m2, ~m3,
                    "Immigration Status", "", "", "",
                    "Naturalized citizen", "-", "-", "-",
                    "Country of Origin", "", "", "",
                    "Ireland", "-", "-", "-",
                    "Education Level", "", "", "",
                    "No high school", "-", "-", "-",
                    "Current Employment", "", "", "",
                    "Currently unemployed", "-", "-", "-",
                    "Total Taxes Paid in 2021", "", "", "",
                    "$500", "-", "-", "-")
  }
  
  attr(rows, "position") <- c(2:3, 6:7, 13:14, 18:19, 23:24)
  
  tab <- modelsummary(mod,
               coef_rename = coef_rename,
               add_rows = rows,
               estimate = "{estimate} ({std.error})",
               statistic = NULL,
               gof_map = c("nobs", "se_type"),
               escape = F,
               output = "markdown")
  
  return(list(fig, tab))
}

# Main results (Figure 2 and Table B.1)
results_study2(bi_pid)

# Results for Independents (Figure E.1)
results_study2(independents)

# Results by education level (Figure E.2)
results_study2(bi_educ)

# Results by household income (Figure E.3)
results_study2(bi_income)

# Results by race (Figure E.4)
results_study2(bi_race)

# Results by nationality (Figure E.5)
results_study2(bi_immstat)

# Results by ideology (Figure E.6)
results_study2(bi_ideo)

# Results by prior immigration attitudes (Figure E.7)
results_study2(bi_immigration)

# Results by Study 1 information treatment (Figure E.8)
results_study2(bi_treatment)

# Results by immigration status in conjoint task (Figure E.9)
results_study2(bi_status, status_subset = TRUE)
              

# Table 3
mod_overall <- lm_robust(selected ~ status + country + skill + job + taxes,
                         clusters = id, weights = weight, data = d2)

pred_overall <- predictions(
  mod_overall, 
  newdata = datagrid(
    status = d2$status, 
    country = d2$country,
    skill = d2$skill,
    job = d2$job,
    taxes = d2$taxes
  )
)

pred_overall %>% 
  tibble() %>%
  filter(status == "Illegal alien"
         & country == "Iraq"
         & skill == "No high school"
         & job == "Currently unemployed" 
         & taxes == "$500") %>%
  select(estimate)

pred_overall %>% 
  tibble() %>%
  filter(status == "Illegal alien"
         & country == "Iraq"
         & skill == "Graduate degree" 
         & job == "Works 60 hours per week" 
         & taxes == "$25,000") %>%
  select(estimate)


pred_overall %>% 
  tibble() %>%
  filter(status == "Naturalized citizen"
         & country == "Ireland" 
         & skill == "No high school"
         & job == "Currently unemployed" 
         & taxes == "$500") %>%
  select(estimate)

pred_overall %>% 
  tibble() %>%
  filter(status == "Naturalized citizen"
         & country == "Ireland" 
         & skill == "Graduate degree" 
         & job == "Works 60 hours per week" 
         & taxes == "$25,000") %>%
  select(estimate)

mod_pid <- lm_robust(selected ~ (status + country + skill + job + taxes) * bi_pid,
                     clusters = id, weights = weight, data = d2)

pred_pid <- predictions(
  mod_pid, 
  newdata = datagrid(
    status = d2$status, 
    country = d2$country,
    skill = d2$skill,
    job = d2$job,
    taxes = d2$taxes,
    bi_pid = d2$bi_pid
  )
)

pred_pid %>% 
  drop_na(bi_pid) %>%
  group_by(bi_pid) %>%
  filter(status == "Illegal alien"
         & country == "Iraq" 
         & skill == "No high school"
         & job == "Currently unemployed" 
         & taxes == "$500") %>%
  ungroup() %>%
  select(bi_pid, estimate)

pred_pid %>% 
  drop_na(bi_pid) %>%
  group_by(bi_pid) %>%
  filter(status == "Illegal alien"
         & country == "Iraq" 
         & skill == "Graduate degree"
         & job == "Works 60 hours per week"
         & taxes == "$25,000") %>%
  ungroup() %>%
  select(bi_pid, estimate)

pred_pid %>% 
  drop_na(bi_pid) %>%
  group_by(bi_pid) %>%
  filter(status == "Naturalized citizen"
         & country == "Ireland" 
         & skill == "No high school"
         & job == "Currently unemployed" 
         & taxes == "$500") %>%
  ungroup() %>%
  select(bi_pid, estimate)

pred_pid %>% 
  drop_na(bi_pid) %>%
  group_by(bi_pid) %>%
  filter(status == "Naturalized citizen"
         & country == "Ireland" 
         & skill == "Graduate degree"
         & job == "Works 60 hours per week"
         & taxes == "$25,000") %>%
  ungroup() %>%
  select(bi_pid, estimate)



# Study 3 results ----------------------------------------------------------------------------------------

results_study3 <- function(subset, xmin = NULL, xmax = NULL, dodge = FALSE) {
  
  mytheme <- theme_tufte(base_size = 12) +
    theme(legend.position = "bottom",
          legend.direction = "horizontal",
          legend.box.margin = margin(-10, -10, -10, -10),
          legend.spacing.y = unit(0, "line"))
  
  p1_overall <- lm_robust(outcome ~ work + assimilation + person, 
                          data = d3, weights = weight, clusters = id) %>%
    tidy() %>% mutate(bi = "Overall")
  
  p1_bi <- d3 %>%
    drop_na({{subset}}) %>%
    group_by({{subset}}) %>%
    group_modify( ~ lm_robust(outcome ~ work + assimilation + person, 
                              data = ., weights = weight, clusters = id) %>%
                    tidy()) %>%
    rename(bi = {{subset}})
  
  if(is.null(xmin)) xmin = -.075
  if(is.null(xmax)) xmax = .175
  
  p1 <- rbind(p1_overall, p1_bi) %>%
    filter(term != "(Intercept)") %>%
    mutate(term = case_when(term == "workHigh" ~ "Economic\ncontributions",
                            term == "personIndividual narrative" ~ "Individual\nnarrative",
                            TRUE ~ "Cultural\nassimilation"),
           term = fct_relevel(term, "Individual\nnarrative"),
           bi = fct_relevel(bi, "Overall")) %>%
    ggplot(aes(estimate, term, 
               color = bi, shape = bi)) +
    geom_vline(xintercept = 0, linetype = "dotted") + 
    geom_point(size = 2, 
               position = position_dodge(-.5)) +
    geom_errorbar(aes(xmin = estimate - 1.96 * std.error,
                      xmax = estimate + 1.96 * std.error),
                  position = position_dodge(-.5),
                  width = 0) +
    geom_errorbar(aes(xmin = estimate - 1.645 * std.error,
                      xmax = estimate + 1.645 * std.error),
                  position = position_dodge(-.5),
                  width = 0, linewidth = 1) +
    labs(y = NULL, color = NULL, shape = NULL,
         x = "Treatment effect estimate with 90% and 95% confidence intervals") +
    scale_x_continuous(limits = c(xmin, xmax), breaks = seq(-.2, .2, .05)) +
    scale_color_grey() +
    mytheme +
    theme(axis.text.y = element_text(hjust=0.5))
  
  p2_overall <- plot_predictions(
    lm_robust(outcome ~ work * assimilation * person, 
              data = d3, weights = weight, clusters = id),
    condition = c("work", "assimilation", "person"), draw = FALSE) %>%
    mutate(bi = "Overall")
  
  p2_bi <- d3 %>%
    drop_na({{subset}}) %>%
    group_by({{subset}}) %>%
    group_modify(
      ~ plot_predictions(
        lm_robust(outcome ~ work * assimilation * person, 
                  data = ., weights = weight, clusters = id),
        condition = c("work", "assimilation", "person"), draw = FALSE)) %>%
    rename(bi = {{subset}})
  
  if(dodge) pos_dodge <- .5
  else pos_dodge <- 0
  
  p2 <- rbind(p2_overall, p2_bi) %>%
    mutate(bi = fct_relevel(bi, "Overall")) %>%
    ggplot(aes(work, estimate,
               color = bi, shape = bi, group = bi)) +
    geom_hline(yintercept = .5, linetype = "dotted") + 
    geom_line(position = position_dodge(pos_dodge)) +  
    geom_point(size = 2, 
               position = position_dodge(pos_dodge)) +
    geom_errorbar(aes(ymin = estimate - 1.96 * std.error,
                      ymax = estimate + 1.96 * std.error),
                  position = position_dodge(pos_dodge),
                  width = 0) +
    geom_errorbar(aes(ymin = estimate - 1.645 * std.error,
                      ymax = estimate + 1.645 * std.error),
                  position = position_dodge(pos_dodge),
                  width = 0, linewidth = 1) +
    scale_color_grey() +
    scale_y_continuous(breaks = seq(.2, .8, .1)) +
    labs(x = "Economic contributions",
         y = "Predicted probability with\n90% and 95% confidence intervals",
         color = NULL, shape = NULL) +
    facet_nested(. ~ person + assimilation) +
    mytheme
  
  fig <- p1 / (p2 + theme(axis.title.y = element_text(margin = margin(r = -100, unit = "pt")))) + 
    plot_annotation(tag_levels = "A") +
    plot_layout(heights = c(1.5, 2))
  
  group1 <- levels(factor(d3 %>% select({{subset}}) %>% pull()))[1]
  group2 <- levels(factor(d3 %>% select({{subset}}) %>% pull()))[2]
  subset1 <- d3 %>% filter({{subset}} == group1)
  subset2 <- d3 %>% filter({{subset}} == group2)

  if (!is.na(group2)) {
    mod_a <- list(
      lm_robust(outcome ~ work + assimilation + person,
                clusters = id, weights = weight, data = d3),
      lm_robust(outcome ~ work + assimilation + person,
                clusters = id, weights = weight, data = subset1),
      lm_robust(outcome ~ work + assimilation + person,
                clusters = id, weights = weight, data = subset2))
    names(mod_a) <- c("Overall", group1, group2)

    mod_b <- list(
      lm_robust(outcome ~ work * assimilation * person,
                clusters = id, weights = weight, data = d3),
      lm_robust(outcome ~ work * assimilation * person,
                clusters = id, weights = weight, data = subset1),
      lm_robust(outcome ~ work * assimilation * person,
                clusters = id, weights = weight, data = subset2))
    names(mod_b) <- c("Overall", group1, group2)
  } else {
    mod_a <- list(
      lm_robust(outcome ~ work + assimilation + person,
                clusters = id, weights = weight, data = d3),
      lm_robust(outcome ~ work + assimilation + person,
                clusters = id, weights = weight, data = subset1))
    names(mod_a) <- c("Overall", group1)

    mod_b <- list(
      lm_robust(outcome ~ work * assimilation * person,
                clusters = id, weights = weight, data = d3),
      lm_robust(outcome ~ work * assimilation * person,
                clusters = id, weights = weight, data = subset1))
    names(mod_b) <- c("Overall", group1)
  }

  tab <- modelsummary(
    list("A" = mod_a, "B" = mod_b),
    coef_map = c("(Intercept)" = "(Intercept)",
                 "workHigh" = "High economic contributions",
                 "assimilationHigh\nassimilation" = "High cultural assimilation",
                 "personIndividual narrative" = "Individual narrative",
                 "workHigh:assimilationHigh\nassimilation" = "High economic × High cultural",
                 "workHigh:personIndividual narrative" = "High economic × Individual",
                 "assimilationHigh\nassimilation:personIndividual narrative" = "High cultural × Individual",
                 "workHigh:assimilationHigh\nassimilation:personIndividual narrative" = "High economic × High cultural × Individual"),
    estimate = "{estimate} ({std.error})",
    statistic = NULL,
    shape = "rcollapse",
    gof_map = c("nobs", "se_type"),
    output = "markdown",
    escape = FALSE)
  
  return(list(fig, tab))
}

# Main results (Figure 3 and Table C.1)
results_study3(bi_pid)

# Results for Independents (Figure C.1)
results_study3(independents, dodge = T, xmin = -.15, xmax = .2)

# Results by education level (Figure C.2)
results_study3(bi_educ, dodge = T, xmin = -.175, xmax = .2)

# Results by household income (Figure C.3)
results_study3(bi_income, dodge = T, xmin = -.13, xmax = .18)

# Results by race (Figure C.4)
results_study3(bi_race, dodge = T, xmin = -.075, xmax = .175)

# Results by ideology (Figure C.5)
results_study3(bi_ideo, xmin = -.125, xmax = .16)


# Table C.2
mod_dems <- lm_robust(outcome ~ work + assimilation + person, weights = weight, clusters = id, data = subset(d3, bi_pid == "Democrats"))
mod_reps <- lm_robust(outcome ~ work + assimilation + person, weights = weight, clusters = id, data = subset(d3, bi_pid == "Republicans"))

clean_terms <- function(tt) {
  tt %>%
    mutate(
      Treatment = case_when(
        str_detect(term, regex("^work", ignore_case = TRUE)) ~ "Economic contributions",
        str_detect(term, regex("^assimilation", ignore_case = TRUE)) ~ "Cultural assimilation",
        str_detect(term, regex("^person", ignore_case = TRUE)) ~ "Individual narrative",
        TRUE ~ NA_character_
      )
    ) %>%
    filter(!is.na(Treatment))
}

td <- tidy(mod_dems) %>% clean_terms() %>%
  select(Treatment, estimate_dem = estimate, se_dem = std.error)

tr <- tidy(mod_reps) %>% clean_terms() %>%
  select(Treatment, estimate_rep = estimate, se_rep = std.error)

results <- td %>%
  left_join(tr, by = "Treatment") %>%
  mutate(
    Democrats   = estimate_dem,
    Republicans = estimate_rep,
    Difference  = Republicans - Democrats,
    SE_diff     = sqrt(se_dem^2 + se_rep^2),               # assumes independent groups
    t           = Difference / SE_diff,
    p           = 2 * (1 - pnorm(abs(t)))                  # two-sided
  ) %>%
  select(Treatment, Democrats, Republicans, Difference, t, p) %>%
  arrange(match(Treatment, c("Economic contributions","Cultural assimilation","Individual narrative")))

results %>% kbl(format = "markdown", digits = 3)




# R version 4.5.2 (2025-10-31 ucrt)
# Platform: x86_64-w64-mingw32/x64
# Running under: Windows 11 x64 (build 26100)
# 
# Matrix products: default
# LAPACK version 3.12.1
# 
# locale:
# [1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8    LC_MONETARY=English_United States.utf8
# [4] LC_NUMERIC=C                           LC_TIME=English_United States.utf8    
# 
# time zone: America/New_York
# tzcode source: internal
# 
# attached base packages:
# [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
# [1] ggh4x_0.3.1            patchwork_1.3.2        ggthemes_5.2.0         kableExtra_1.4.0       marginaleffects_0.31.0
# [6] modelsummary_2.5.0     broom.helpers_1.22.0   estimatr_1.0.6         lubridate_1.9.4        forcats_1.0.1         
# [11] stringr_1.6.0          dplyr_1.1.4            purrr_1.2.0            readr_2.1.5            tidyr_1.3.1           
# [16] tibble_3.3.0           ggplot2_4.0.1          tidyverse_2.0.0       
# 
# loaded via a namespace (and not attached):
# [1] gtable_0.3.6        bayestestR_0.17.0   xfun_0.54           insight_1.4.4       tzdb_0.5.0          vctrs_0.6.5        
# [7] tools_4.5.2         generics_0.1.4      datawizard_1.3.0    parallel_4.5.2      pkgconfig_2.0.3     tinytable_0.15.1   
# [13] data.table_1.17.8   checkmate_2.3.3     RColorBrewer_1.1-3  S7_0.2.0            lifecycle_1.0.4     compiler_4.5.2     
# [19] farver_2.1.2        textshaping_1.0.4   codetools_0.2-20    htmltools_0.5.8.1   yaml_2.3.10         Formula_1.2-5      
# [25] pillar_1.11.1       parallelly_1.46.0   tidyselect_1.2.1    digest_0.6.38       performance_0.15.3  future_1.68.0      
# [31] stringi_1.8.7       listenv_0.10.0      labeling_0.4.3      labelled_2.16.0     fastmap_1.2.0       grid_4.5.2         
# [37] cli_3.6.5           magrittr_2.0.4      cards_0.7.1         utf8_1.2.6          future.apply_1.20.1 broom_1.0.11       
# [43] withr_3.0.2         scales_1.4.0        backports_1.5.0     timechange_0.3.0    rmarkdown_2.30      httr_1.4.7         
# [49] globals_0.18.0      hms_1.1.4           evaluate_1.0.5      knitr_1.50          haven_2.5.5         parameters_0.28.3  
# [55] viridisLite_0.4.2   rlang_1.1.6         Rcpp_1.1.0          glue_1.8.0          xml2_1.5.1          effectsize_1.0.1   
# [61] svglite_2.2.2       rstudioapi_0.17.1   R6_2.6.1            tables_0.9.33       systemfonts_1.3.1   texreg_1.39.4    